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We show that a Gaussian Process model can be combined with a small number (of order 100) of 
scattering calculations to provide a multi-dimensional dependence of scattering observables on the 
experimentally controllable parameters (such as the collision energy or temperature) as well as the 
potential energy surface (PES) parameters. For the case of Ar - CeHg collisions, we show that 
200 classical trajectory calculations are sufficient to provide a 10-dimensional hypersurface, giving 
the dependence of the collision lifetimes on the collision energy, internal temperature and 8 PES 
parameters. This can be used for solving the inverse scattering problem, the efficient calculation 
of thermally averaged observables, for reducing the error of the molecular dynamics calculations by 
averaging over the PES variations, and the analysis of the sensitivity of the observables to individual 
parameters determining the PES. Trained by a combination of classical and quantum calculations, 
the model provides an accurate description of the quantum scattering cross sections, even near 
scattering resonances. 


The reliable scattering calculations of dynamical prop¬ 
erties of molecules are required in almost any research 
field related to molecular physics. In particular, the ex¬ 
periments on collisional cooling of molecules to cold and 
ultracold temperatures [1]], chemical reaction dynamics 
0, the development of new pressure standards 0, 0], 
astrophysics and astrochemistry [jj rely on accurate cal¬ 
culations of molecular collision cross sections. Currently, 
there are two major problems with the ab initio calcula¬ 
tions of molecular dynamics observables. The first prob¬ 
lem is the inaccuracy of the potential energy surfaces 
(PES). Unfortunately, even the most sophisticated quan¬ 
tum chemistry calculations produce the PES with un¬ 
certainties that lead to significant (and often unknown) 
errors in the dynamical calculations. This sensitivity to 
PES inaccuracies is especially detrimental for low tem¬ 
perature applications (cold molecules, ultracold chem¬ 
istry, astrophysics and pressure standards) 0-Q. The 
second problem is related to the numerical complexity of 
the quantum dynamics calculations 00- For complex 
molecules with many degrees of freedom, accurate dy¬ 
namical calculations are extremely time-consuming and 
it is often impossible to compute enough results for ac¬ 
curate averaging over the collision or internal energies of 
the colliding partners. 

In the present work we propose a solution to these two 
problems. In order to account for the PES uncertainties, 
the dynamical results can be averaged over variations of 
the PES. If the computed observables are averaged over 
variations of each individual PES parameter, producing 
an expectation interval of the observables, the ab initio 
dynamical calculations can have fully predictive power 
(with error bars). However, the outcome of a molecu¬ 
lar collision is generally a complicated function of many 
(ten or more) PES parameters. It is impossible to ob¬ 
tain the dependence of the collision observables on the 
individual PES parameters by the direct scattering cal¬ 
culations. We show that such a dependence can be ob¬ 


tained by a combination of a small number (on the order 
of 100) of scattering calculations with a Gaussian Process 
(GP) model llj,[l2|. We show that the same model can 
be used to obtain the accurate dependence of the scat¬ 
tering observables on the collision or internal energies of 
the molecules, with a small number of scattering calcu¬ 
lations. The result is an accurate global dependence of 
the scattering observables on the collision energy, inter¬ 
nal energy and every individual parameter of the PES 
surface. This global dependence can be used to average 
the computed observables over variations of the individ¬ 
ual PES parameters, as well as over the collision and 
internal energies in order to produce thermally average 
observables. It can also be used to analyze the influence 
of the individual PES parameters on the scattering out¬ 
come. This makes the model proposed here a unique tool 
for the analysis of the effects of the PES topology on the 
molecular scattering dynamics. 


Widely used in engineering technologies 00 , the 
GP model can be viewed as a technique for interpolation 
in a multi-dimensional space. We choose the GP model 
because it is an efficient non-parametric method. There 
is no need to fit data by analytical functions so the model 
is expected to work for any distribution of scattering ob¬ 
servables and to become more accurate when trained by 
more computed observables. Given the scattering ob¬ 
servables computed at a small number of randomly cho¬ 
sen points in the multi-dimensional parameter space, the 
GP model learns from correlations between the values 
of these scattering observables to produce a smooth de¬ 
pendence on all the underlying parameters. As an illus¬ 
trative example, we consider the scattering of benzene 
molecules CeH 6 by rare gas (Rg) atoms He - Xe. The 
PES surface for CeH 6 - Rg interactions is characterized 
by 8 parameters. We consider two scattering observables 
151-1171]: the collision lifetimes and the scattering cross 
sections. We address the following questions: how many 
scattering calculations are sufficient to train a GP model 
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to produce an accurate global dependence on all the un¬ 
derlying parameters? Can the GP model be used to make 
predictions of the scattering observables for one collision 
system based on the known properties of another colli¬ 
sion system? Can the GP model be used to characterize 
the scattering observables near quantum resonances? 

We consider a scattering observable O as a function 
of q parameters described by vector x. The components 
of the vector x = (xi,X 2 , ■■■ ,x q ) T can be the collision 
energy, the internal energy and/or the parameters rep¬ 
resenting the PES. We assume that O is known from a 
classical or quantum dynamics computation at a small 
number of x values. Our first goal is to construct an ef¬ 
ficient model that, given a finite set of 0(x), produces a 
global dependence of the scattering observable on s. If 
the observable is known from a measurement or a rigor¬ 
ous quantum calculation as a function of some param¬ 
eters Xi - e.g., the collision energy - we show that the 
model can be adjusted to produce the global dependence 
of O on x that reproduces the accurate data, even if the 
dynamical calculation method is inaccurate. 

We assume that the scattering observable of interest at 
any a; is a realization of a Gaussian process F(-), char¬ 
acterized by a mean function /i(-), constant variance cr 2 
and correlation function R(-, •). For any fixed x , F(x) 
is a value of a function randomly drawn from a family 
of functions Gaussian-distributed around /r(-). Conse¬ 
quently, the multiple outputs F(x) and F(x') at x and 
x' jointly follow a multivariate normal distribution de¬ 
fined by /x(-), cr 2 , and i?(-, •) 0, EH- We assume the 
following form for the correlation function I20I - 23 1: 


R(x, x') = exp < — 


w»Fi - Xu 


'ip 


( 1 ) 


and write 


F(x) = ^2 hj(x)/3j + Z(x) = h (x) T /3 + Z{x), (2) 
i=i 


where h = ( h 1 (x ),..., hk(x)) is a vector of k regression 
functions [24|, (3 = (/3i,/? 2,-'' ,/?fc) T is a vector of un¬ 
known coefficients, and Z{-) is a Gaussian random func¬ 
tion with zero mean. The problem is thus reduced to 
finding (3, p and L2 = (wi, u> 2 , ■ ■ • , w g ) T . 

We spread n input vectors xi,...,x n evenly through¬ 
out a region of interest and compute the desired ob¬ 
servable O at each x,; with a classical or quantum dy¬ 
namics method. The outputs of a GP at these points 


Y n = 


^F(xi), F(x 2 ), ■ ■ ■ ,F(xn2j follow a multivari¬ 


ate normal distribution with the mean vector H/3 and 
the covariance matrix er 2 A. Here, H is a n x k de¬ 
sign matrix with ith row filled with the k regressors 
hi(xi), h, 2 (xi), ■ ■ ■ ,hk(xi ) at site aq , and Aisanxn 
matrix with the elements A (i,j) = R(xi, xj). 


Given fi, the maximum likelihood estimators (MLE) 
of (3 and cr 2 have closed-form solutions m 

/3(H) = (H t A" 1 H)- 1 H t A^ 1 Y n (3) 


a 2 (n) =-<Y n -U/3) T A~ 1 (Y n -U/3) (4) 

n 

To find the MLE of fi, we fix p and maximize the log- 
likelihood function 

log£(n|l™) = -i [n\oga 2 + log(det(A)) + n] (5) 

numerically by an iterative computation of the determi¬ 
nant | A| and the matrix inverse A -1 . 

The goal is to make a prediction of the scattering ob¬ 
servable at an arbitrary x = Xq. Because the values 
To = F(x 0 ) at a; 0 and the outputs at training sites are 
jointly distributed, the conditional distribution of possi¬ 
ble values To = F{x 0 ) given the values Y n is a normal 
distribution with the conditional mean and variance 


m(x 0 )* = h{x 0 ) T f3 + AjA" 1 ( Y n - H/3) (6) 

a* 2 (x 0 ) = ( r 2 (l-A ( [A^ 1 Ao), (7) 


where A 0 = (R(x 0l aq), R(x 0 , x 2 ), ■■■, R(x 0: x n )) T ^ is 
specified by the now known correlation function 
Eq. © provides the GP model prediction for the value 
of the scattering observable at Xq 

To illustrate the applicability and accuracy of the GP 
model, we first compute the collision lifetimes of benzene 
molecules with Rg atoms [US]. We use the classical 
trajectory (CT) method described in Ref. (25j. As shown 
in Ref. [28|, the CeH 6 - Rg PES can be expressed as a 
sum over terms describing the interaction of Rg with the 
C-C and C-H bond fragments, characterized by 8 pa¬ 
rameters. We first fix the PES parameters to describe 
the C 6 H 6 - Ar system and focus on the dependence of 
the lifetimes on two parameters: the collision energy E 
and the rotational temperature T r . Figure 1 shows the 
results of the CT calculations illustrating that the colli¬ 
sion lifetime exhibits an inverse correlation with E, while 
no apparent correlation with T r . Figure 1 (c) shows the 
global surface of the lifetime as a function of E and T r 
obtained from the GP model with hi = 1, /ii>i = 0 and p 
set to 1.95. To quantify the prediction accuracy of the GP 

model, we calculate the errors £e = \J ^ ~ Vi) 2 


and £s = £e/ (f/max - 2/min), where y t are the computed 
values and yt are the GP model predictions. For the 
model with only 20 scattering calculations used as train¬ 
ing points, £e = 9.36 ps and £s = 7.93 %. With the 
number of the scattering calculations increased to 50, the 
errors decrease to £e = 5.17 ps and £s = 4.38 %. 

The scattering calculations presented in the upper pan¬ 
els of Figure 1 cannot be interpreted to assume any sim¬ 
ple functional form. In addition, the vastly different gra¬ 
dients of the T r and E dependence may make the conclu¬ 
sions based on calculations at fixed values of one of the 
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FIG. 1. (a) and (b): The lifetime dependence on the rota¬ 

tional temperature and collision energy for CeHg - Ar colli¬ 
sions. (c): The surface produced by the GP model. The lines 
connect the values (circles) computed from the classical tra¬ 
jectories with the values predicted by the GP model. (d): The 
surface produced by the GP model for CgHg - Rg collision 
lifetimes vs the atomic mass and the PES depth for T r = 4 K 
and E = 4 cm -1 . The surface (c) is produced with only 20 
scattering calculations on input and has the normalized error 
£s < 8 %. The surface (d) is produced with 40 scattering 
calculations and has the error es = 5.09 %. 


parameters misleading. In contrast, the surface plot in 
Figure 1(c) clearly illustrates that the collision lifetimes 
decrease monotonically with both T r and E. The effect 
of the rotational temperature is much weaker especially 
when E > 5 cm -1 and there is no strong two-way inter¬ 
action between T r and E. The GP model surface can be 
used to evaluate thermally averaged collision lifetimes by 
integrating the FI-dependence at given T r . 

The GP model can be extended to multiple collision 
systems for the predictions of the collision properties of 
a specific collision system based on the known collision 
properties of another system. To illustrate this, we con¬ 
sider the lifetimes of the long-lived complexes formed by 
benzene in collisions with Rg atoms He - Xe. As the col¬ 
lision system is changed from CeH 6 - He to CeH 6 - Xe, 
there are two varying factors that determine the change 
of the collision dynamics: the reduced mass and the PES. 

As before, we use the GP model F(x) = /3 + Z(x), 
with x now representing the atomic mass ha and the in¬ 
teraction strength D e at the global minimum of the atom 
- molecule PES obtained by scaling the Ar - CeH 6 PES. 


We fix T r = 4 K and E = 4 cm -1 , and compute the colli¬ 
sion lifetimes at 40 randomly chosen points in the interval 
of ha and D e [4g/mol, 130g/mol] x [80cm _1 ,520cm _1 ], 
which covers all of the Rg - CgHg systems. These 40 
calculation points are then used to train the GP model 
to produce the surface plot shown in Figure 1 (d). The 
error £5 of the surface is 5.09 %. The plot reveals that in¬ 
creasing both ha and D e enhances the collision lifetimes 
and that the reduced-mass dependence of the collision 
lifetimes is very weak compared to the dependence on 
the interaction strength. 



Computed Lifetime (ps) 


FIG. 2. Accuracy of the GP model with variable PES param¬ 
eters for the prediction of the collision lifetimes. The scatter 
plot compares the predicted values with the computed values. 
The error of the GP model is the deviation of the points from 
the diagonal line. This GP model is trained by only 200 
scattering calculations, enough to produce a 10-dimensional 
hypersurface with the error Es = 4 %. Left inset: Energy de¬ 
pendence of the collision lifetime for Ar - CgH6 with the error 
interval obtained by varying all the individual PES parame¬ 
ters by ±3 %. Right inset: Relative effect of the variation of 
T r , E and the PES parameters on the collision lifetimes. The 
filled area of the bars shows the uncorrelated contribution of 
the corresponding variable and the open area - the effect that 
depends on one or more other variables. 

The GP model can be exploited to explore the role of 
the individual PES parameters on the observables. To 
illustrate this, we now consider that x contains 8 param¬ 
eters giving the analytical form of the Rg - CeH 6 PES 
[ 2 ^ . in addition to E and T r . We calculate the lifetimes 
at 200 randomly selected points in this parameter space 
and use these points to train the GP model. Figure 2 
compares the predicted values with the calculated values 
for another set of 70 randomly selected points. The plot 
corresponds to the model error £s = 4 %. 

The 10-parameter GP model contains a wealth of in- 
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formation on the dependence O(x). For example, one 
can perform a sensitivity analysis by using the functional 
analysis of variance decomposition 29j-31| to determine, 
which of the PES parameters have the strongest impact 
on the observable (right inset of Figure 2). Of the 8 PES 
parameters, the location of the potential well due to the 
interactions of Rg with the C-C bonds for the parallel 
approach [28| is the most important factor determining 
the collision lifetime. The model can also use be used 
to compute the uncertainties due to global variation of 
the PES. Figure 2 (left inset) shows the interval of the 
lifetimes obtained by the simultaneous ±3 % variation of 
all 8 PES parameters. 

We now consider the applicability of the GP model 
to quantum scattering calculations. The quantum re¬ 
sults are often affected by resonances |2[, [32J, leading to 
wild variations of the scattering observables in a small 
range of the underlying parameters. If applied directly 
to such the case, the GP model is unstable because steep 
variation of the correlations leads to singularities in A” 1 
[33|. This is illustrated in Figure 3, showing the GP 
model predictions trained directly by 60 quantum cal¬ 
culations of cross sections for rotationally inelastic He - 
CgHg scattering, randomly chosen at E between 1 and 10 
cm -1 . The instability of the GP model arises from the 
wild variations of the scattering cross sections near reso¬ 
nances. We repeated these calculations for the elastic and 
state-resolved rotationally inelastic cross sections shown 
in Figure 4 (a-c) of Ref. |27j. In each case, we found 
that the wild variation of the quantum results leads to 
unstable GP model predictions. 

However, the GP model can be extended to model the 
time-consuming quantum scattering calculations with 
the help of efficient classical dynamics calculations. To 
do this, we introduce a more complex GP as (mJ 


E(-) = pF(.) + G(-)+e, 


( 8 ) 


where F(-) and G(-) are independent Gaussian random 
functions, with G(-) characterizing the difference between 
the CT and QM calculations and effectively describing 
the inaccuracy of the classical trajectory method. The 
calculations are performed in two steps. First, the CT 
calculations are used to train the GP model F(-). In 
the second step, the QM and CT calculations are used 
together to train the model G(-) in Eq. ©, using the 
parameters of F(-) and treating p and e as variable pa¬ 
rameters. This fixes the models F(-) and G(-) as well as 
p and e. 

The accuracy of this combined quantum - classical 
model is illustrated in Figure 3, showing that the model 
provides an accurate energy dependence of the cross sec¬ 
tions, even near scattering resonances. The CT calcu¬ 
lations in a two-function model © stabilize the model, 
removing the errors arising from the resonant variation 
of the quantum results. We applied the two-step model 
© to the calculations for the elastic and state-resolved 


rotationally inelastic cross sections shown in Figure 4 (a- 



FIG. 3. GP models (solid curves) of quantum scattering 
cross sections (symbols) for CeH6 - He collisions. Blue dashed 
line: quantum calculations are used directly to train the GP 
model ©. Red solid line: A combination of classical and 
quantum results is used in a hybrid GP model lf8)l . The CT 
results stabilize the GP model predictions of the quantum cal¬ 
culations. The models are trained by the points represented 
by squares. The circles are used to illustrate the accuracy. 

In summary, we have shown that a Gaussian Pro¬ 
cess model combined with a small number of scattering 
calculations can be used to obtain an accurate multi¬ 
dimensional dependence of the scattering observables on 
the experimentally controllable parameters and the PES 
parameters. Specifically, we showed that the GP model 
trained only by 20 CT calculations produces a depen¬ 
dence of the CeHg - Ar collision lifetimes on the colli¬ 
sion energy and the rotational temperature of benzene, 
with the normalized error eg < 8 %. Trained by 200 
calculations, the GP model produces a 10-dimensional 
dependence of the collision lifetimes on the collision en¬ 
ergy, the rotational temperature and 8 individual PES 
prameters, with the error £s < 4 %. We have introduced 
a hybrid GP model that can be trained by a combina¬ 
tion of classical and quantum dynamics calculations in 
order to model the quantum results. We showed that 
this model works even in the vicnity of quantum scat¬ 
tering resonances, where the direct fit of the quantum 
results by means of a GP model is unstable. The models 
described here are expected to find a wide range of appli¬ 
cations, from fitting the interaction potentials by solving 
the inverse scattering problem, to analyzing the depen¬ 
dence of scattering observations on external parameters, 
to calibrating the accuracy of the scattering calculation 
methods. For example, the inverse scattering problem 
can be approached with the help of Eq. ©, where F(-) 
is parametrized by unknown PES parameters and E(-) 
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models the experimental data. The best estimates of 
the unknown PES parameters can then be found by a 
Markov-chain Monte Carlo method [35| _, in a procedure 
similar to one recently applied in Ref. [361 ]. 

We thank Dr. Zhiying Li for allowing us to use her 
codes for the classical and quantum dynamics calcula¬ 
tions. This work is supported by NSERC of Canada. 
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